library(ggplot2)
library(vegan)
library(dplyr)
library(purrr)
library(readr)
library(stringr)
library(plotly)
## A set of colors
color.list <- c("#E64B35FF", "#4DBBD5FF", "#00A087FF", "#3C5488FF", "#F39B7FFF", "#8491B4FF",
"#91D1C2FF", "#B09C85FF", "#FAFD7CFF", "#82491EFF", "#B7E4F9FF", "#FB6467FF",
"#526E2DFF", "#E762D7FF", "#FAE48BFF", "#A6EEE6FF", "#95CC5EFF")
## ggplot2 preset theme
figtheme <- theme_classic() +
theme(text = element_text(size=23,face='bold'),
axis.title.y=element_text(margin=margin(0,15,0,0)),axis.title.x=element_text(margin=margin(15,0,0,0)),
plot.margin = unit(c(1,1,1,1), "cm"),
plot.title = element_text(margin=margin(0,0,15,0), hjust=0.5))
theme_set(figtheme)
## Auxilary function to merge two profiles
merge2 <- function(x, y){
aux <- function(x) read_tsv(x, col_names=c("tax", str_remove(basename(x), '[\\._\\-].*')))
if(is.character(x)){
dat.x <- aux(x)
}else{
dat.x <- x
}
dat.y <- aux(y)
full_join(dat.x, dat.y)
}
profile.list <- read_lines(params$profile_list)
tax.profile <- reduce(profile.list, merge2)
tax.profile <- tibble::column_to_rownames(tax.profile, 'tax')
tax.profile[is.na(tax.profile)|tax.profile<params$presence_absence_thre] <- 0 ## impute NAs and set low abundance taxon to 0
tax.profile <- tax.profile[rowSums(tax.profile) != 0, ] ## remove empty taxa
tax.profile
distance <- vegdist(t(tax.profile))
cmds <- cmdscale(distance, eig = TRUE)
perc <- (cmds$eig/sum(cmds$eig))[1:2]*100
plot.dat <- data.frame(cmds$points)
p <- ggplot(plot.dat, aes(x=X1, y=X2)) +
geom_density_2d(color='grey') +
geom_point(size=3) +
scale_color_manual(values=color.list[c(1,4)]) +
labs(x=sprintf('PCoA 1 [%.1f%%]', perc[1]), y=sprintf('PCoA 2 [%.1f%%]', perc[2]))
ggplotly(p)
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: CentOS release 6.5 (Final)
##
## Matrix products: default
## BLAS/LAPACK: /mnt/software/unstowable/miniconda3-4.6.14/envs/shotgunMetagenomics_r_v3.6.0/lib/libopenblasp-r0.3.7.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] plotly_4.9.1 stringr_1.4.0 readr_1.3.1 purrr_0.3.2
## [5] dplyr_0.8.0.1 vegan_2.5-6 lattice_0.20-38 permute_0.9-5
## [9] ggplot2_3.1.1
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.5 xfun_0.11 splines_3.6.1 colorspace_1.4-1
## [5] htmltools_0.3.6 viridisLite_0.3.0 yaml_2.2.0 mgcv_1.8-28
## [9] rlang_0.3.4 pillar_1.3.1 later_1.0.0 glue_1.3.1
## [13] withr_2.1.2 plyr_1.8.4 munsell_0.5.0 gtable_0.3.0
## [17] htmlwidgets_1.5.1 evaluate_0.13 labeling_0.3 knitr_1.26
## [21] httpuv_1.5.2 crosstalk_1.0.0 parallel_3.6.1 Rcpp_1.0.1
## [25] xtable_1.8-4 scales_1.0.0 promises_1.1.0 jsonlite_1.6
## [29] mime_0.6 hms_0.4.2 digest_0.6.18 stringi_1.4.3
## [33] shiny_1.3.2 grid_3.6.1 tools_3.6.1 magrittr_1.5
## [37] lazyeval_0.2.2 tibble_2.1.1 cluster_2.0.8 crayon_1.3.4
## [41] tidyr_0.8.3 pkgconfig_2.0.2 MASS_7.3-51.3 Matrix_1.2-17
## [45] data.table_1.12.6 assertthat_0.2.1 rmarkdown_1.12 httr_1.4.0
## [49] R6_2.4.0 nlme_3.1-139 compiler_3.6.1